Changes in the gradient percolation transition caused by an Allee effect 



Michael T. Gastncr, 1,2 Beata Oborny, 3 Alcxey B. Ryabov, 1 and Bcrnd Blasius 1 

institute for the Chemistry and Biology of the Marine Environment, 
Carl von Ossietzky Universitat, Carl-von-Ossietzky-Strafie 9-11, 26111 Oldenburg, Germany 

2 Department of Mathematics, Complexity and Networks Progamme, 
Imperial College London, South Kensington Campus, London SW7 2AZ, United Kingdom 
Department of Plant Taxonomy and Ecology, Lordnd Eotvos University, 
Pdzmdny Peter stny. 1/C, H-1117, Budapest, Hungary 

The establishment and spreading of biological populations depends crucially on population growth 
at low densities. The Allee effect is a problem in those populations where the per-capita growth rate 
at low densities is reduced. We examine stochastic spatial models in which the reproduction rate 
changes across a gradient g so that the population undergoes a 2D-percolation transition. Without 
the Allee effect, the transition is continuous and the width w of the hull scales as in conventional (i.e., 



uncorrelated) gradient percolation, 
is first order and w oc g~ 2(i . 



w oc g 



However, with a strong Allee effect the transition 
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It is not just human relationships which obey the 
rule "two's company, three's a crowd." Negative den- 
sity dependence, defined as a decrease of the per-capita 
growth rate with increasing population density, is com- 
mon among almost all species at high densities, where 
overcrowding and the depletion of resources limit further 
growth. The most common model for negative density 
dependence is the logistic equation which assumes that 
the per-capita growth rate decreases linearly with the 
population size P, 
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where t is time, r is the intrinsic rate of increase, and K 
the carrying capacity. 

If r, K > 0, Eq.[T]is characterized by a negative density 
dependence for all population sizes P. For some small 
populations, however, a positive density dependence can 
be observed. The latter is called a demographic Allee 
effect, named after Warder Clyde Allee, who described it 
first and supported the theory with examples from var- 
ious animal species from insects to mammals Small 
populations can suffer from reduced growth rates for 
various reasons. Frequently, a collective behavior (e.g., 
defense against predators) becomes inefficient when the 
group is small. Additionally, small populations are less 
efficient in modifying the environment to their own ben- 
efit. For example, plant individuals in aggregations can 
reduce frost or desiccation, but only when the density in 
the clump is sufficiently high 0]. To generalize density 
dependence, Volterra proposed to replace the right-hand 
side of Eq. [T] with a quadratic function of P 0] , 



1 dP 
P~dt 



= -A + BP-CP 2 , A 7 B 7 C>0. 



(2) 



If B 2 — 4AC > 0, Eq.[2]has two stable fixed points, unlike 
Eq. [TJ which has only one, so that the long-term behav- 
ior of Eq. [2] depends on the initial population density. If 



P> (B- \/B 2 - AAC)/2C at t = 0, the population will 
approach a positive limit, whereas a smaller initial popu- 
lation will become extinct. Several other formulations of 
the Allee effect have been suggested in the past decades, 
some including stochastic and spatial effects US (see 
chapter 3.5 in Ref. @ for a review). They all have in 
common that a strongly positive density dependence ac- 
celerates the extinction of small populations. 

The work described here is motivated by the question: 
what are the consequences of an Alice effect on popula- 
tions that live at a margin of a geographic range? Be- 
cause such populations usually have low densities, one 
can expect that it matters greatly for the success of es- 
tablishment and spreading if an Allee effect is present 
or not d, 0]. In this Letter, we investigate the situa- 
tion near a geographic margin with two models where 
the density changes across space from low to high values. 
We show that a strong Alice effect makes the percolation 
transition at the margin discontinuous, causing scaling 
behavior different from previously studied types of gra- 
dient percolation [1, Q . 

Our models are stochastic cellular automata whose lo- 
cal rules correspond to discretized versions of Eq. Q] or 
[3] Both cellular automata operate on a two-dimensional 
honeycomb lattice where the sites are either populated 
(A) or vacant (0, Fig. QJ,). They can change their state 
by local death and birth events. In both models, deaths 
are Poisson processes: 

• A — > 0: A populated site becomes vacant with rate 
1. 

In our first model, the rate, with which a vacant site 
becomes populated by a local birth event, is exactly pro- 
portional to the number of neighbors: 

• A — > 2A: A vacant site at position (x, y) with k 
populated adjacent sites becomes itself populated 
at the rate b(x) • fc/6. 
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FIG. 1: (a) Sites in the spatial models are placed in the 
centers of the hexagons in a honeycomb lattice. Gray cells 
represent populated (A), white cells vacant (0) sites. In this 
example, the focal site in the center has k = 4 populated 
neighbors. If this site is populated, it will die during a small 
time interval dt with probability equal to dt. If, on the other 
hand, this site is vacant, it will become populated with prob- 
ability (fc/6) b(x f )dt in the GCP or [\k(k - 1)/15] b(x f )dt in 
the GAP, where b(xf) is the local birth attempt rate, (b) 
Typical snapshot of the GCP, (c) of the GAP. Dark gray: 
the largest populated cluster. Light gray: all other populated 
sites. Black curve: percolation hull. The mean hull position 
x and the width of the fluctuations w are indicated at the 
bottom. 



The second model implements a local Allee effect by re- 
quiring at least one pair of neighbors for successful births. 
The rule A — > 2A is replaced with: 

• 2A — > 3A: A vacant site at (x, y) with k neighbors 
[i.e., |fc(fc — 1) pairs of neighbors] becomes popu- 
lated with rate b(x) ■ ^k(k — 1)/15. 
The denominators 6 and 15 are the maximum number of 
neighbors and the maximum number of neighbor pairs, 
respectively. Sites are updated in a random order with 
the rates stated above following the algorithm of Ref. [IcJ 
The function b(x) can be interpreted as the rate with 
which an individual in column x attempts to produce 
offspring on an adjacent site. A birth attempt succeeds 
only if that site is vacant. In the case of 2A — > 3A, 
success further depends on a second neighbor adjacent to 
the newly born individual. If b(x) is a constant, then the 
first model is equivalent to a contact process [llj , and our 
second model becomes a special case of "Schldgl's second 
model" Ell. 



Our work differs from previous studies by assuming 
a constant gradient g > in the birth attempt rate, 
b(x) = gx. Long-range gradients are important in ecol- 
ogy because the environmental conditions can change 
gradually over distances larger than the distance of dis- 
persal within one generation (e.g., along a hillside or 
across geographic latitudes). We call A — > 2 A a gra- 
dient contact process (GCP) [13, E3 and 2A -> 3A a 
gradient Allee process (GAP). As long as the initial pop- 
ulation density is sufficiently high in regions with high 
birth rates, the population density reaches a steady state 
independent of the initial conditions (see supplement). 

Because g > 0, the steady-state density of populated 
sites grows in both the GCP and the GAP as x, and 
hence b, increases. At small x, the populated sites form 
small isolated patches (light gray sites in Figs. [TJd and 
[TJ;) whereas at large x most populated sites belong to 
one large cluster (dark gray). The curve along which the 
largest cluster touches the largest contiguous vacant area 
(black curve) is the percolation hull @ [3]. If the pop- 
ulated sites provide habitat or food for another species, 
the hull marks the borderline between the connected and 
fragmented occurrence of this resource. An example is 
a treeline across an altitudinal or latitudinal gradient. 
Births and deaths cause the position and shape of the 
percolation hull to fluctuate. The average position of the 
hull x and the characteristic width of the fluctuations w 
depend on g and the model (GCP versus GAP). We com- 
pute x and w as the mean and the standard deviation of 
the distribution of x-coordinates along the hull during 
several independent runs. 

In Fig. [lib), the number of sites in the GCP's largest 
cluster increases gradually from left to right. The in- 
crease is much more abrupt in the GAP (Fig. [TJ;) which 
generates fewer isolated clusters. This impression can be 
confirmed by looking at local densities in the transition 
region near x. Because w is the relevant length scale in 
this region, we investigate subsystems located between 
x — w and x + w (Fig. [2ji) . We determine the cluster 
that has the largest number of sites N within a 2w x 2w 
square. In both models, the mean of N scales approxi- 
mately as w Df with Df = 91/48, the fractal dimension 
of the incipient infinite cluster in standard (i.e., uncorre- 
lated, grad ient-free, nondirectcd) two-dimensional perco- 
lation [15J (see supplement). But the distribution of the 
N sites is different in the two models. 

To see this quantitatively, let us define the cluster den- 
sity p as N divided by the number of all (populated or 
vacant) sites in the square. The distributions of p, aggre- 
gated over independent runs of the GCP and the GAP, at 
different y-coordinates and at different times, are shown 
in Fig. [2Jd. The GCP distribution has a single peak at 
intermediate densities whereas the GAP has two local 
maxima, one at zero and another one at high density. In 
analogy to thermodynamics, where a bimodal probability 
distribution of an order parameter is an indication of a 
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FIG. 2: (a) To distinguish between a continuous and a first- 
order transition, we investigate squares of dimension 2w x 2w 
centered at x. We measure the fraction of sites p that belong 
to the cluster covering the largest area within the square. 
(Typically, though not necessarily, this is the largest cluster 
on the entire lattice, shown in dark gray), (b) From data 
of several independent runs, we obtain the probability distri- 
bution for p, represented as a histogram. We show data for 
g — 5 ■ 10~ 5 . The smaller the gradient, the more weight is 
concentrated in the two peaks of the GAP distribution, (c) 
The distribution of cluster sizes s in the stripe \x—x\ < w. We 
include clusters which are at least partially in this stripe, but 
exclude the system's largest cluster. The dashed line oc s _T 
is the tangent to the GCP distribution. The inset shows the 
data collapse for the GCP. The exponents are those of stan- 
dard 2D-percolation: r=^,a=|f,i/=| []]|. 



first-order phase transition [16[ , percolation in the GAP 
can be interpreted as a first-order transition between two 
steady states of either zero or of a positive density. Thus 
only a population larger than a critical density is able 
to grow and, when this density is exceeded, the cluster 
size grows abruptly, reminiscent of recent reports of "ex- 
plosive percolation" [13], but generated here by a purely 
local rule. 

In some explosive percolation models, cluster sizes were 
recently shown to follow power-law distributions, casting 
doubts on whether the transition is truly first order . 
In the GAP, by contrast, the cluster size distribution 
Pcs(s) in the stripe \x — x\ < w gives further evidence 
in support of a first-order transition (Fig. [2t). While the 
GCP distribution follows the scaling behavior expected 
for two-dimensional percolation with a continuous transi- 
tion p cs (s) = s" T / cs (sg 1/[ ' T(l/+1)I ), the GAP distribution 
does not show indications of scaling. There is neither a 
power-law decay nor a dependence on the gradient even 
for large cluster sizes. Instead, consistent with a first- 
order transition, there is a characteristic size for the GAP 
clusters. 
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FIG. 3: The hull width to as a function of the gradient g. 
The lines are least-squares fits to the data. Error bars are 
smaller than the symbol sizes. 



Another percolation property affected by the Allee ef- 
fect is the scaling relation between w and g (Fig. [3]). For 
the GCP, we find w oc .g- aGCP , a GC p = 0.572(3) (95% 
confidence interval) . The GAP width also follows a power 
law, but with a smaller exponent oqap = 0.26(1). (The 
same exponents are observed if the hull is replaced with 
the accessible external perimeter (l9| . see supplement.) 
The exponent oqcp is in agreement with vj (y+1) = 4/7, 
which was derived for uncorrelated percolation based on 
the divergence of the correlation length 0. There is 
no analogous relation for ogap, because the correlation 
length is finite at a first-order transition. The result 
that w, nevertheless, scales with g in the GAP is sur- 
prising, considering that scaling in stochastic gradient 
models has so far always been linked to divergent corre- 
lation lengths 0] . 

Although the spatial width of the hull increases with 
decreasing g, the transition zone becomes, in terms of 
the birth rate b, more confined. This becomes clear by 
plotting pi c {b), the probability that a site with birth rate 
b belongs to the largest cluster (Fig. g}. I n both models, 
Pic(b) approaches a limiting function as g — > + with a 
sharp increase at the percolation thresholds 6 p ,gcp = 
2.260(1) and 6 PiG ap = 7.7340(3). For finite g, p lc (b) 
obeys the finite-size scaling pi c (b,g) = g c fic(\b — b p \g a ~ 1 ), 
where a is the hull width exponent (accp or ogap) and 
fi c is a modcl-dcpcndcnt scaling function. In the GCP, 
we expect c = j3/(u + 1) = 5/84. For the GAP, however, 
the first-order transition demands that pi c has a discon- 
tinuity in the limit g — > + and therefore cqap = 0. The 
insets of Fig. 2] show a remarkable data collapse for the 
anticipated exponents. 

Why is percolation in the GAP unconventional? Let 
us denote by n(b, t) the probability that a site with birth 
attempt rate b is populated at time t. The mean-field 
equations for n to lowest order in g are 

GCP: ^=- n + b(l-n)n + ^b(l-n)^, (3) 
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FIG. 4: The probability p; c that a site belongs to the largest 
cluster as a function of the site's birth attempt rate b in (a) 
the GCP, (b) the GAP for various gradients g. Insets: The 
functions collapse if the coordinates are rescaled. In both in- 
sets, the same nine data sets are shown as in the main panels. 



GAP: ^ = -n + b(l -n)n 2 + 
at 



10 



b(l-n) 



5n 



db 2 



On 
db 



(4) 



where b(x) = gx (see supplement for the derivation and 
numerical solutions). If g = 0, Eq. [3] turns into the lo- 
gistic equation ((T|) with r = b — 1, P = bKn/r, and 
Eq. H becomes Eq. [2] with P = n, A = 1, B = C = b. 
For g — > + , Eq. [3] has a continuous stationary solution: 
n = max(0, 1 — l/b). In Eq. |4] on the other hand, n 
develops a discontinuity in the limit g — > + where it 
suddenly jumps from zero ton~ 0.77. If the percolation 
threshold is at a probability n p — 0.5, as in uncorrclatcd 
honeycomb lattices [15j . the GCP hull is at a position 
where n(b) varies smoothly. The GAP hull, by contrast, 
lies at a position where n(b) changes abruptly. 

In summary, the GCP and the GAP behave fundamen- 
tally differently near the margin of the populated range. 
The GCP, a model without any Allee effect, possesses 
the same characteristic features as previously reported 
for uncorrelated gradient percolation [Toj | . The Allee ef- 



fect in the GAP changes the situation drastically: the 
percolation transition is first order and the hull width 
diverges more slowly for g — > + . 
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